##==============================================================##
## Title: Leader Motivated Behavior: Voting by Mail in the      ##
##         2020 General Election                                ##
## Authors: Seth C. McKee, Enrijeta Shino, and Daniel A. Smith  ##
## Journal: Public Opinion Quarterly                            ##
##==============================================================##

rm(list = ls())

# List of required packages
packages <- c("dplyr", "tidyr", "lme4", "lubridate", "formattable", 
              "scales", "stringr", "lemon", "janitor", "readxl", 
              "readr", "sf", "tigris", "ggplot2", "sandwich", "broom",
              "lmtest", "purrr", "texreg", "rstudioapi", 
              "descr", "forcats", "scales", "lmtest", "sandwich",
              "ggeffects", "emmeans", "multcomp", "gridExtra", "texreg",
              "readstata13", "data.table")

# Install missing packages
install_if_missing <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
  }
}

# Apply the installation function to all packages
invisible(lapply(packages, install_if_missing))

# Load the packages
lapply(packages, library, character.only = TRUE)


#=================#
# CES 2016 SURVEY #
#=================#

load("/Users/Enrijeta/Dropbox/POQ_VBM_CES/ces_data/POQ_Replication/ces16.RData")


#====================#
## RECODE VARIABLES ##
#====================#

# Self-reported: Voted in 2016 election: 1 = yes; 0 = no
ces16$voted2016 <- ifelse(ces16$CC16_401 == "I definitely voted in the General Election.", 1, 0)
table(ces16$voted2016)
prop.table(table(ces16$voted2016))


# Validated vote 2016: 1= voted/unknown method; 0 = missing
ces16$CL_E2016GVM1 <- ces16$CL_E2016GVM
ces16$CL_E2016GVM1[ces16$CL_E2016GVM1 == ""] <- 99
ces16$vote16_validated <- ifelse(ces16$CL_E2016GVM1 == 99, 0, 1)
table(ces16$vote16_validated)
prop.table(table(ces16$vote16_validated))


# Validated vote method: 0 = in-person (early + ED), 1 = mail/absentee, unknown dropped
ces16$votemethod_validated <- with(ces16, ifelse(CL_E2016GVM== "earlyVote" | CL_E2016GVM == "polling", 0, 
                                          ifelse(CL_E2016GVM == "absentee" | CL_E2016GVM == "mail", 1, NA)))

table(ces16$votemethod_validated)
prop.table(table(ces16$votemethod_validated))


# Self-reported vote method: 0 = in-person (early + ED), 1 = mail
ces16$votemethod_r <- with(ces16, ifelse((CC16_403 == "In person on election day" |
                                          CC16_403 == "In person before election day (early)"), 0,
                                  ifelse(CC16_403 == "Voted by mail (or absentee)", 1, NA)))

table(ces16$votemethod_r)
prop.table(table(ces16$votemethod_r))


# Self-reported vote method for validated voters only: 0 = in-person (early + ED), 1 = mail
ces16$votemethod_rv <- with(ces16, ifelse((CC16_403 == "In person on election day" | 
                                           CC16_403 == "In person before election day (early)") & vote16_validated == 1, 0, 
                                   ifelse(CC16_403 == "Voted by mail (or absentee)" & vote16_validated == 1, 1, NA)))

table(ces16$votemethod_rv)
prop.table(table(ces16$votemethod_rv))


# VBM voting misreport, validated: 1 = misreport; 0 = correctly reported VBM
ces16$vbm_misreport <- with(ces16, ifelse(votemethod_validated == 1 & votemethod_rv != 1, 1,
                                   ifelse(votemethod_validated == 1 & votemethod_rv == 1, 0, NA)))
table(ces16$vbm_misreport)


# Voted in 2016, self-reported and validated measure: 1 = validated voter; 0 = otherwise
ces16$voter2016 <- ifelse(ces16$voted2016 == 1 & ces16$vote16_validated== 1, 1, 0)
table(ces16$voter2016, useNA = "ifany")


# Education: less than hs; high school edu; some college; college grad; post grad)
ces16$edulevel <- with(ces16, ifelse(educ == "No HS", 1,
                              ifelse(educ == "High school graduate", 2,
                              ifelse(educ == "Some college" | educ == "2-year", 3, 
                              ifelse(educ == "4-year", 4,
                              ifelse(educ == "Post-grad", 5, NA))))))
table(ces16$edulevel)


# Race: 0 = white, 1 = black, 2 = Hispanic, 3 = race other 
ces16$race_cat <- with(ces16, ifelse(race == "White", 0,
                              ifelse(race == "Black", 1,
                              ifelse(race == "Hispanic", 2, 3))))
table(ces16$race_cat)
ces16$race_cat <- factor(ces16$race_cat)


# Gender: 0 = male; 1 = female
ces16$female <- ifelse(ces16$gender == "Female", 1, 0)
table(ces16$female)


# Income continuous
ces16$faminc[ces16$faminc == "Prefer not to say"] <- NA

ces16$income <- as.numeric(ces16$faminc)
table(ces16$income)


# Social staus: 1 = married; 0 = not married
ces16$married <- ifelse(ces16$marstat == "Married", 1, 0)
table(ces16$married)


# Age in 2016
ces16$age <- 2016 - ces16$birthyr
table(ces16$age)


# News interest: 1 = no interest --> 4 = a great deal
ces16$news <- with(ces16, ifelse(newsint == "Hardly at all", 1,
                          ifelse(newsint == "Only now and then", 2, 
                          ifelse(newsint == "Some of the time", 3,
                          ifelse(newsint == "Most of the time", 4, NA)))))
table(ces16$news)


# Vote choice 2016: 1 = Trump, 0 = Other
ces16$trump <- with(ces16, ifelse(CC16_364c == "Donald Trump (Republican)", 1, 0))
table(ces16$trump)


############################################################################################
# 2016: 28 states
# NCSL's NO excuse for absentee ballot states: AK, AZ, CA, DC, FL, GA, HI, ID, IL, 
#  IA, KS, ME, MD, MN, MT, NE, NV, NJ, NM, NC, ND, OH, OK, SD, UT, VT, WI, WY
# Source: http://web.archive.org/web/20161011024719/https:/www.ncsl.org/research/elections-and-campaigns/absentee-and-early-voting.aspx
############################################################################################

ces16$noexcuse_vbm <- with(ces16, ifelse(inputstate=="Alaska" | inputstate=="Arizona" | inputstate=="California" | 
                                         inputstate=="District of Columbia" | inputstate=="Florida" | inputstate=="Georgia" |
                                         inputstate=="Hawaii" | inputstate=="Idaho" | inputstate=="Illinois" | inputstate=="Iowa" | 
                                         inputstate=="Kansas" | inputstate=="Maine" | inputstate=="Maryland" | inputstate=="Minnesota" | 
                                         inputstate=="Montana"| inputstate=="Nebraska" | inputstate=="Nevada" | inputstate=="New Jersey" | 
                                         inputstate=="New Mexico" | inputstate=="North Carolina" | inputstate=="North Dakota" | 
                                         inputstate=="Ohio" | inputstate=="Oklahoma" | inputstate=="South Dakota" | inputstate=="Utah" |
                                         inputstate=="Vermont" | inputstate=="Wisconsin" | inputstate=="Wyoming", 1, 0)) 

table(ces16$noexcuse_vbm)


# Dummies for all mail election states: CO, OR, WA
# https://www.ncsl.org/research/elections-and-campaigns/vopp-table-18-states-with-all-mail-elections.aspx
ces16$allmail <- ifelse(ces16$inputstate == "Colorado" | ces16$inputstate == "Oregon"| ces16$inputstate == "Washington", 1, 0)
table(ces16$allmail)


# Battleground states: 1= yes; 0 = no 
# Colorado, Florida, Iowa, Maine, Michigan, Nebraska, Nevada, New Hampshire, North Carolina, 
# Ohio, Pennsylvania, Virginia, Wisconsin (13 states)
ces16$battleground <- with(ces16, ifelse(inputstate=="Colorado"| inputstate=="Florida"| inputstate=="Iowa"| inputstate=="Maine"| 
                                         inputstate=="Michigan"| (inputstate=="Nebraska" & cdid115 == "2")| inputstate=="Nevada"| 
                                         inputstate=="New Hampshire"| inputstate=="North Carolina"| inputstate=="Ohio"| 
                                         inputstate=="Pennsylvania"| inputstate=="Virginia"| inputstate=="Wisconsin", 1, 0))
table(ces16$battleground)


# Party ID: 1 = Ind/other/not sure; 2 = Republican; 3 = Democrat
ces16$pid_3cat <- with(ces16, ifelse(pid3 == "Democrat", 3, 
                              ifelse(pid3 == "Republican", 2, 1)))

table(ces16$pid_3cat)
ces16$pid_3cat <- factor(ces16$pid_3cat)


# ideology: 1= moderate/not sure, 2= conservative, 3= liberal
ces16$ideo3 <- with(ces16, ifelse(ideo5 == "Liberal" | ideo5 == "Very liberal", 3, 
                           ifelse(ideo5 == "Conservative" | ideo5 == "Very conservative", 2, 1)))

table(ces16$ideo3)
ces16$ideo3 <- factor(ces16$ideo3)


# Church attendance: 1 never ---> 6 more than once a week
ces16$pew_churatd1 <- as.numeric(ces16$pew_churatd)
ces16$pew_churatd1[ces16$pew_churatd1 == 7] <- NA
ces16$church <- 7 - ces16$pew_churatd1
table(ces16$church)


# year
ces16$year <- 0


##====================##
# Descriptives Table 1 #
##====================##

ces16 %>%
  filter(!is.na(votemethod_validated), !is.na(trump)) %>%
  tabyl(votemethod_validated, trump, weights = commonweight_vv) %>%
  adorn_totals("row") %>%
  adorn_percentages("col") %>%
  adorn_pct_formatting()


ces16 %>%
  filter(!is.na(votemethod_rv), !is.na(trump)) %>%
  tabyl(votemethod_rv, trump, weights = commonweight_vv) %>%
  adorn_totals("row") %>%
  adorn_percentages("col") %>%
  adorn_pct_formatting()


# Correlation coefficient for Trump support and Republican
ces16$pid_3cat1 <- as.numeric(ces16$pid_3cat)
ces16$pid_3cat1[ces16$pid_3cat1==3 | ces16$pid_3cat1==1 ] <- 0
cor(ces16$trump, ces16$pid_3cat1, use = "na.or.complete")  # 0.543


##=============================##
# OLS Regression Models Table 2 #
##=============================##

# model 1
m2016_v <- lm(votemethod_validated ~ trump + pid_3cat + ideo3 + female + edulevel + 
                                     race_cat + age + married + income + church + news + 
                                     noexcuse_vbm + allmail + battleground, data = ces16)
summary(m2016_v)

# model 2
m2016_s <- lm(votemethod_rv ~ trump + pid_3cat + ideo3 + female + edulevel + 
                              race_cat + age + married + income + church + news + 
                              noexcuse_vbm + allmail + battleground, data = ces16)
summary(m2016_s)

# model 3
m2016 <- lm(vbm_misreport ~ trump + pid_3cat + ideo3 + female + edulevel + 
                            race_cat + age + married + income + church + news + 
                            noexcuse_vbm + allmail + battleground, data = ces16)
summary(m2016)


#=====================================================#
# Marginal Effects for Vote Method Misreport, Plot 1a #
#=====================================================#
ces16 <- data.frame(ces16)

# Validated vote method marginal effects
validated_trump16 <- ggpredict(m2016_v, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 1, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

validated_trump16$model <- "Validated"


# Self-reported vote method marginal effects
selfreport_trump16 <- ggpredict(m2016_s, terms = c("trump"), 
                                typical = "mode", condition = c(pid_3cat = 1, female = 0, race_cat = 0, 
                                                                married = 1, allmail = 0, battleground = 1), 
                                ci.lvl = 0.95, vcov.fun = "vcovCL")

selfreport_trump16$model <- "Self-reported"


# Misreported vote method marginal effects
misreport_trump16 <- ggpredict(m2016, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 1, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

misreport_trump16$model <- "Misreported"

ces16 <- data.table(ces16)


## Combine data 
plot_data16 <- rbind(validated_trump16, selfreport_trump16, misreport_trump16)


#=========#
# Plot 1a #
#=========#

# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.5) # move them .05 to the left and right

# Make them factors to preserve order in x-axis as desired. 
plot_data16$model <- factor(plot_data16$model, levels = c("Validated","Self-reported","Misreported"))


plot2016 <- ggplot(plot_data16, aes(x = x, y = predicted, colour = model)) + 
  geom_point(aes(shape = model), position = pd, size = 3.5) + 
  guides(shape = FALSE) + 
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high, colour = model, width = 0.0),
    size = 1, position = pd
  ) + 
  geom_hline(yintercept = 0.2372990, linetype = "dashed", color = "red", size = 0.8) +
  scale_y_continuous(
    breaks = seq(0, 0.55, by = 0.1),
    limits = c(0, 0.55),
    expand = c(0, 0)
  ) + 
  scale_x_continuous(
    breaks = c(0, 1),
    labels = c("0" = "Non Trump\nSupporter", "1" = "Trump\nSupporter"),
    expand = c(0.05, 0)
  ) +
  xlab("") +
  ylab("Predicted values") +
  scale_shape_manual(values = c(15, 17, 18)) +
  guides(colour = guide_legend(override.aes = list(shape = c(15, 17, 18)))) +
  theme_bw() +
  theme(
    plot.margin = unit(c(0.7, 0.7, 0.7, 0.7), "cm"),
    text = element_text(size = 18),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    axis.title = element_text(size = 16),
    axis.text.y = element_text(size = 15),
    axis.text.x = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 12)
  ) +
  scale_color_discrete(name = " ")

# Make it black and white
plot2016 + scale_colour_manual(
  name = " ",
  values = c("Validated" = "#000000", "Self-reported" = "#000000", "Misreported" = "#000000")
)


# save graph on overleaf
ggsave(
  filename = "/Users/Enrijeta/Dropbox/Apps/Overleaf/Leader Motivated Behavior/plots/plot16.eps",
  width = 5,
  height = 6,
  device = "eps"
)


##=================##
##=================##
## CES 2020 SURVEY ##
##=================##
##=================##

load("/Users/Enrijeta/Dropbox/POQ_VBM_CES/ces_data/POQ_Replication/ces20.RData")

#====================#
## RECODE VARIABLES ##
#====================#

# Self-reported: Voted in 2020 election: 1 = yes; 0 = no
ces20$voted2020 <- ifelse(ces20$CC20_401 == 5, 1, 0)
table(ces20$voted2020, useNA = "ifany")
prop.table(table(ces20$voted2020))


# Validated vote 2020: 1= voted/unknown method; 0 = missing
ces20$CL_2020gvm1 <- ces20$CL_2020gvm
ces20$CL_2020gvm1[is.na(ces20$CL_2020gvm1)] <- 99
ces20$vote20_validated <- ifelse(ces20$CL_2020gvm1 == 99, 0, 1)
table(ces20$vote20_validated)
prop.table(table(ces20$vote20_validated))


# Validated vote method: 0 = in-person (early + ED), 1 = mail/absentee, (unknown, dropped)
ces20$votemethod_validated <- with(ces20, ifelse(CL_2020gvm== 1 | CL_2020gvm == 3, 1, 
                                          ifelse(CL_2020gvm == 2 | CL_2020gvm == 4, 0, NA)))

table(ces20$votemethod_validated)
prop.table(table(ces20$votemethod_validated))


# Self-reported vote method: 0 = in-person (early + ED), 1 = mail; (idk answers dropped)
ces20$votemethod_r <- with(ces20, ifelse((CC20_403 == 1 | CC20_403 == 2), 0, 
                                  ifelse(CC20_403 == 3, 1, NA)))

table(ces20$votemethod_r)
prop.table(table(ces20$votemethod_r))


# Self-reported vote method, validated voters: 0 = in-person (early + ED), 1 = mail; (idk answers dropped)
ces20$votemethod_rv <- with(ces20, ifelse((CC20_403 == 1 | CC20_403 == 2) & vote20_validated == 1, 0, 
                                   ifelse(CC20_403 == 3 & vote20_validated == 1, 1, NA)))

table(ces20$votemethod_rv)
prop.table(table(ces20$votemethod_rv))


# VBM voting misreport, validated: 1 = misreport; 0 = correctly reported VBM
ces20$vbm_misreport <- with(ces20, ifelse(votemethod_validated == 1 & votemethod_rv != 1, 1,
                                   ifelse(votemethod_validated == 1 & votemethod_rv == 1, 0, NA)))
table(ces20$vbm_misreport)


# Voted in 2020, self-reported and validated measure: 1 = validated voter; 0 = otherwise
ces20$voter2020 <- ifelse(ces20$voted2020 == 1 & ces20$vote20_validated== 1, 1, 0)
table(ces20$voter2020, useNA = "ifany")


# Voted in 2020, misreported turnout (self-reported and validated measure): 1 =  misreported overerport; 0 = otherwise
ces20$voter2020_misreport <- ifelse(ces20$voted2020 == 1 & ces20$vote20_validated== 0, 1, 0)
table(ces20$voter2020_misreport, useNA = "ifany")


# Education: less than hs; high school edu; some college; college grad; post grad)
ces20$edulevel <- with(ces20, ifelse(educ == 1, 1,
                              ifelse(educ == 2, 2,
                              ifelse(educ == 3 | educ == 4, 3, 
                              ifelse(educ == 5, 4,
                              ifelse(educ == 6, 5, NA))))))
table(ces20$edulevel)


# Race: 0 = white, 1 = black, 2 = Hispanic, 3 = race other 
ces20$race_cat <- with(ces20, ifelse(race == 1, 0,
                              ifelse(race == 2, 1,
                              ifelse(race == 3, 2, 3))))

table(ces20$race_cat)
ces20$race_cat <- factor(ces20$race_cat)


# Gender: 0 = male; 1 = female
ces20$female <- ifelse(ces20$gender == 2, 1, 0)
table(ces20$female)


# Income continuous
ces20$faminc_new[ces20$faminc_new == 97] <- NA

ces20$income <- ces20$faminc_new
table(ces20$income)


# 1: married; 0: not married
ces20$married <- ifelse(ces20$marstat == 1, 1, 0)
table(ces20$married)


# Age in 2020
ces20$age <- 2020 - ces20$birthyr
table(ces20$age)


# Vote choice 2020: 1 = Trump, 0 = Other candidate
ces20$trump <- with(ces20, ifelse(CC20_364b == 1, 1, 0))
table(ces20$trump)


# News interest: 1= no interest --> 4 a great deal
ces20$newsint[ces20$newsint == 7] <- NA
ces20$news <- 5 - ces20$newsint
table(ces20$news)


# COVID-19
ces20$covid19 <- with(ces20, ifelse(CC20_309a_1 == 1 | CC20_309a_2 == 1 | 
                                    CC20_309a_3 == 1| CC20_309a_4 == 1, 1,
                             ifelse(CC20_309a_5 == 1, 0, 0)))

table(ces20$covid19)


###########################################################################################################
# 2020: 35 states (Added temporary no-excuse DE, MA, SC, WV)
# NCSL's NO excuse abs ballot: AK, AZ, CA, DC, DE, FL, GA,  ID, IL, IA, KS, MA,
#  ME, MD, MI, MN, MT, NE, NV, NJ, NM, NC, ND, OH, OK, PA, RI, SD, SC, UT, VT, VA,
#  WA, WV, WI, WY; https://www.ncsl.org/research/elections-and-campaigns/absentee-and-early-voting.aspx
###########################################################################################################

ces20$noexcuse_vbm <- with(ces20, ifelse(inputstate==2 | inputstate==4 | inputstate==6 | inputstate==10 | inputstate==11 | 
                                         inputstate==12 | inputstate==13 | inputstate==25 | inputstate==16 | 
                                         inputstate==17 | inputstate==19 | inputstate==20 | inputstate==23 | 
                                         inputstate==24 | inputstate==26 | inputstate==27 | inputstate==30 | 
                                         inputstate==31 | inputstate==32 | inputstate==34 | inputstate==35 |
                                         inputstate==37 | inputstate==38 | inputstate==39 | inputstate==40 | 
                                         inputstate==54 | inputstate==42 | inputstate==44 | inputstate==46 | 
                                         inputstate==49 | inputstate==50 | inputstate==51 | inputstate==45 | 
                                         inputstate==55 | inputstate==56, 1, 0))

table(ces20$noexcuse_vbm)


# Dummies for all mail election states: CO, OR, HI, WA
# https://www.ncsl.org/research/elections-and-campaigns/vopp-table-18-states-with-all-mail-elections.aspx

ces20$allmail <- ifelse(ces20$inputstate == 8 | ces20$inputstate == 15| 
                        ces20$inputstate == 53| ces20$inputstate == 41, 1, 0)
table(ces20$allmail)


# Battleground states: 1= yes; 0 = no 
# Arizona, Florida, Georgia, Maine CD2, Michigan, Minnesota, Nebraska CD2, Nevada, 
# New Hampshire, North Carolina, Pennsylvania, Wisconsin (12 states)
ces20$battleground <- with(ces20, ifelse(inputstate==4| inputstate==12| inputstate==13| (inputstate==23 & cdid116 == "2")| 
                                         inputstate==26| inputstate==27| (inputstate==31 & cdid116 == "2") | 
                                         inputstate==32| inputstate==33| inputstate==37| inputstate==42| 
                                         inputstate==55, 1, 0))
table(ces20$battleground)


# Party ID: 1 = Ind/other/not sure; 2 = Republican; 3 = Democrat
ces20$pid_3cat <- with(ces20, ifelse(pid3 == 1, 3, 
                              ifelse(pid3 == 2, 2, 1)))

table(ces20$pid_3cat)
ces20$pid_3cat <- factor(ces20$pid_3cat)


# Ideology: 1= moderate/not sure, 2= conservative, 3= liberal
ces20$ideo3 <- with(ces20, ifelse(ideo5 == 1 | ideo5 == 2, 3, 
                           ifelse(ideo5 == 4 | ideo5 == 5, 2, 1)))

table(ces20$ideo3)
ces20$ideo3 <- factor(ces20$ideo3)


# Church attendance: 1 never ---> 6 more than once a week
ces20$pew_churatd[ces20$pew_churatd == 7] <- NA
ces20$church <- 7 - ces20$pew_churatd
table(ces20$church)


# Year
ces20$year <- 1



##====================##
# Descriptives Table 1 #
##====================##

ces20 %>%
  filter(!is.na(votemethod_validated), !is.na(trump)) %>%
  tabyl(votemethod_validated, trump, weights = vvweight) %>%
  adorn_totals("row") %>%
  adorn_percentages("col") %>%
  adorn_pct_formatting()


ces20 %>%
  filter(!is.na(votemethod_rv), !is.na(trump)) %>%
  tabyl(votemethod_rv, trump, weights = vvweight) %>%
  adorn_totals("row") %>%
  adorn_percentages("col") %>%
  adorn_pct_formatting()


# Correlation coefficient for Trump support and Republican
ces20$pid_3cat1 <- as.numeric(ces20$pid_3cat)
ces20$pid_3cat1[ces20$pid_3cat1==3 | ces20$pid_3cat1==1 ] <- 0
cor(ces20$trump, ces20$pid_3cat1, use = "na.or.complete") # 0.631



##=======================================================##
# OLS Regression Models for Vote Method Misreport Table 2 #
##=======================================================##

m2020_v <- lm(votemethod_validated ~ trump + pid_3cat + ideo3 + female + edulevel + 
                                     race_cat + age + married + income + church + news + 
                                     noexcuse_vbm + allmail + battleground, data = ces20)
summary(m2020_v)


m2020_s <- lm(votemethod_rv ~ trump + pid_3cat + ideo3 + female + edulevel + 
                              race_cat + age + married + income + church + news + 
                              noexcuse_vbm + allmail + battleground, data = ces20)
summary(m2020_s)



m2020 <- lm(vbm_misreport ~ trump + pid_3cat + ideo3 + female + edulevel + 
                            race_cat + age + married + income + church + news + 
                            noexcuse_vbm + allmail + battleground, data = ces20)
summary(m2020)



#====================================================#
# Marginal Effects for Vote Method Misreport Plot 1b #
#====================================================#
ces20 <- data.frame(ces20)

# Validated vote method marginal effects
validated_trump20 <- ggpredict(m2020_v, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 1, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

validated_trump20$model <- "Validated"


# Self-reported vote method marginal effects
selfreport_trump20 <- ggpredict(m2020_s, terms = c("trump"), 
                                typical = "mode", condition = c(pid_3cat = 1, female = 0, race_cat = 0, 
                                                                married = 1, allmail = 0, battleground = 1), 
                                ci.lvl = 0.95, vcov.fun = "vcovCL")

selfreport_trump20$model <- "Self-reported"


# Misreported vote method marginal effects
misreport_trump20 <- ggpredict(m2020, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 1, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

misreport_trump20$model <- "Misreported"

ces20 <- data.table(ces20)


### Combine data 

plot_data20 <- rbind(validated_trump20, selfreport_trump20, misreport_trump20)


#=========#
# Plot 1b #
#=========#

# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.5) # move them .05 to the left and right

plot_data20$model <- factor(plot_data20$model, levels = c("Validated","Self-reported","Misreported"))


plot2020 <- ggplot(plot_data20, aes(x = x, y = predicted, colour = model)) + 
  geom_point(aes(shape = model), position = pd, size = 3.5) + 
  guides(shape = FALSE) + 
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high, colour = model, width = 0.0),
    size = 1, position = pd
  ) + 
  geom_hline(yintercept = 0.2372990, linetype = "dashed", color = "red", size = 0.8) +
  scale_y_continuous(
    breaks = seq(0, 0.55, by = 0.1),
    limits = c(0, 0.55),
    expand = c(0, 0)
  ) + 
  scale_x_continuous(
    breaks = c(0, 1),
    labels = c("0" = "Non Trump\nSupporter", "1" = "Trump\nSupporter"),
    expand = c(0.05, 0)
  ) +
  xlab("") +
  ylab("Predicted values") +
  scale_shape_manual(values = c(15, 17, 18)) +
  guides(colour = guide_legend(override.aes = list(shape = c(15, 17, 18)))) +
  theme_bw() +
  theme(
    plot.margin = unit(c(0.7, 0.7, 0.7, 0.7), "cm"),
    text = element_text(size = 18),
    panel.grid.minor = element_blank(),
    legend.position = "bottom",
    axis.title = element_text(size = 16),
    axis.text.y = element_text(size = 15),
    axis.text.x = element_text(size = 15),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 12)
  ) +
  scale_color_discrete(name = " ")

# Make it black and white
plot2020 + scale_colour_manual(
  name = " ",
  values = c("Validated" = "#000000", "Self-reported" = "#000000", "Misreported" = "#000000")
)


# save graph on overleaf
ggsave(
  filename = "/Users/Enrijeta/Dropbox/Apps/Overleaf/Leader Motivated Behavior/plots/plot20.eps",
  width = 5,
  height = 6,
  device = "eps"
)


#=========================================#
#=========================================#
# CES 2016 and 2020 Pooled Models Table 2 #
#=========================================#
#=========================================#

# Convert ces16 to a data.table 
ces16 <- as.data.table(ces16)
ces20 <- as.data.table(ces20)


# Select specific columns
ces161 <- ces16[, .(votemethod_validated, votemethod_rv, votemethod_r, vbm_misreport, trump, 
                    pid_3cat, ideo3, female, edulevel, race_cat, age, married, 
                    income, church, news, noexcuse_vbm, allmail, battleground, year, inputstate)]


ces201 <- ces20[, .(votemethod_validated, votemethod_rv, votemethod_r, vbm_misreport, trump, pid_3cat, ideo3, 
                    female, edulevel, race_cat, age, married, income, church, news, noexcuse_vbm,
                    allmail, battleground, year, inputstate)] 


# combine
ces_all <- rbind(ces161, ces201)


##=============================##
# OLS Regression Models Table 2 #
##=============================##

all_v <- lm(votemethod_validated ~  trump + pid_3cat + ideo3 + female + edulevel + 
                                    race_cat + age + married + income + church + news + 
                                    noexcuse_vbm + allmail + battleground + year, data = ces_all)
summary(all_v)


all_s <- lm(votemethod_rv ~  trump + pid_3cat + ideo3 + female + edulevel + 
                             race_cat + age + married + income + church + news + 
                             noexcuse_vbm + allmail + battleground + year, data = ces_all)
summary(all_s)


all <- lm(vbm_misreport ~ trump + pid_3cat + ideo3 + female + edulevel + 
                          race_cat + age + married + income + church + news + 
                          noexcuse_vbm + allmail + battleground + year, data = ces_all)
summary(all)



##=============================##
# DiD Regression Models Table 2 #
##=============================##

all_vd <- lm(votemethod_validated ~  trump*year + pid_3cat + ideo3 + female + edulevel + 
                                     race_cat + age + married + income + church + news + 
                                     noexcuse_vbm + allmail + battleground, data = ces_all)
summary(all_vd)


all_sd <- lm(votemethod_rv ~  trump*year + pid_3cat + ideo3 + female + edulevel + 
                              race_cat + age + married + income + church + news + 
                              noexcuse_vbm + allmail + battleground, data = ces_all)
summary(all_sd)


alld <- lm(vbm_misreport ~ trump*year + pid_3cat + ideo3 + female + edulevel + 
                           race_cat + age + married + income + church + news + 
                           noexcuse_vbm + allmail + battleground, data = ces_all)
summary(alld)


#================#
# Table 2 Output #
#================#

# variable names
map <- list("trump" = "Trump supporter",
            "pid_3cat2" = "Republican",
            "pid_3cat3" = "Democrat",
            "noexcuse_vbm" = "No excuse absentee",
            "allmail" = "All mail",
            "battleground" = "Battleground",
            "year" = "Year 2020",
            "trump:year" = "TrumpXYear 2020",
            "ideo32" = "Conservative",
            "ideo33" = "Liberal",
            "female" = "Female",
            "edulevel" = "Education",
            "race_cat1" = "Black",
            "race_cat2" = "Hispanic",
            "race_cat3" = "Other",
            "age" = "Age",
            "married" = "Married",
            "income" = "Income",
            "church" = "Church attendance",
            "news" = "Political interest",
            "(Intercept)" = "(Intercept)")

# Specify list of models 
models = list(m2016_v, m2020_v, all_v, all_vd,
              m2016_s, m2020_s, all_s, all_sd,
              m2016,m2020, all, alld)

# Specify caption 
caption = "Estimating VBM in 2016 and 2020: Validated VBM, Self-reported VBM, and Misreported VBM"


texreg(models,
       custom.coef.map = map,
       caption = caption,
       caption.above = TRUE,
       digits = 3,
       dcolumn = FALSE,
       stars = numeric(0), 
       override.se = list(
         summary(m2016_v)$coef[,4],
         summary(m2020_v)$coef[,4],
         summary(all_v)$coef[,4],
         summary(all_vd)$coef[,4],
         summary(m2016_s)$coef[,4],
         summary(m2020_s)$coef[,4],
         summary(all_s)$coef[,4],
         summary(all_sd)$coef[,4],
         summary(m2016)$coef[,4],
         summary(m2020)$coef[,4],
         summary(all)$coef[,4],
         summary(alld)$coef[,4]
       )
)




#======================================================================#
# Appendix: Marginal Effects for Vote Method Misreport for Republicans #
#======================================================================#
ces16 <- data.frame(ces16)

# Validated vote method marginal effects
validated_trump16 <- ggpredict(m2016_v, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 2, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

validated_trump16$model <- "Validated"


# Self-reported vote method marginal effects
selfreport_trump16 <- ggpredict(m2016_s, terms = c("trump"), 
                                typical = "mode", condition = c(pid_3cat = 2, female = 0, race_cat = 0, 
                                                                married = 1, allmail = 0, battleground = 1), 
                                ci.lvl = 0.95, vcov.fun = "vcovCL")

selfreport_trump16$model <- "Self-reported"



# Misreported vote method marginal effects
misreport_trump16 <- ggpredict(m2016, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 2, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

misreport_trump16$model <- "Misreported"

ces16 <- data.table(ces16)


### Combine data 

plot_data16 <- rbind(validated_trump16, selfreport_trump16, misreport_trump16)


# Plot 1

# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.5) # move them .05 to the left and right

# Make them factors to preserve order in x-axis as desired. 
plot_data16$model <- factor(plot_data16$model, levels = c("Validated","Self-reported","Misreported"))


plot2016 <- ggplot(plot_data16, aes(x = x, y = predicted, colour = model)) + 
  geom_point(aes(shape = model), position = pd, size=3.5) + 
  guides(shape=FALSE) + 
  geom_errorbar(
    aes(
      ymin = conf.low, ymax = conf.high,
      colour = model, width = 0.0
    ), size=1, position = pd) + theme_light() + 
  geom_hline(yintercept = 0.2170260, linetype="dashed",
             color = "red", size = 0.8) +
  scale_y_continuous(breaks=seq(0, 0.55, by = 0.1), limits=c(0, 0.55)) + 
  scale_x_continuous(breaks=c(0, 1), 
                     labels = c("0" = "Non Trump Supporter", "1" = "Trump Supporter"), expand = c(0, 0))  +
  xlab("") +
  ylab("Predicted values for Republicans") + 
  scale_shape_manual(values = c(15, 17, 18)) + 
  guides(colour = guide_legend(override.aes = list(shape = c(15, 17, 18)))) +
  theme_bw() + theme(text = element_text(size = 15), 
                     panel.grid.minor = element_blank(), legend.position = "bottom", 
                     axis.title = element_text(size = 12),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12),
                     legend.text = element_text(size = 12),
                     legend.title=element_text(size=10)) +
  scale_color_discrete(name=" ") 

plot2016

# make it black and white 
plot2016 + scale_colour_manual(name=" ", 
                               values = c("Validated" = "#000000", "Self-reported" = "#000000", "Misreported" = "#000000"))


# save graph on overleaf
ggsave(filename = "/Users/Enrijeta/Dropbox/Apps/Overleaf/Leader Motivated Behavior/plots/plot16_rep.pdf", width = 6, height = 6)


#======================================================#
# Appendix: Marginal Effects for Vote Method Misreport #
#======================================================#
ces20 <- data.frame(ces20)

# Validated vote method marginal effects
validated_trump20 <- ggpredict(m2020_v, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 2, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

validated_trump20$model <- "Validated"


# Self-reported vote method marginal effects
selfreport_trump20 <- ggpredict(m2020_s, terms = c("trump"), 
                                typical = "mode", condition = c(pid_3cat = 2, female = 0, race_cat = 0, 
                                                                married = 1, allmail = 0, battleground = 1), 
                                ci.lvl = 0.95, vcov.fun = "vcovCL")

selfreport_trump20$model <- "Self-reported"


# Misreported vote method marginal effects
misreport_trump20 <- ggpredict(m2020, terms = c("trump"), 
                               typical = "mode", condition = c(pid_3cat = 2, female = 0, race_cat = 0, 
                                                               married = 1, allmail = 0, battleground = 1), 
                               ci.lvl = 0.95, vcov.fun = "vcovCL")

misreport_trump20$model <- "Misreported"

ces20 <- data.table(ces20)


### Combine data 

plot_data20 <- rbind(validated_trump20, selfreport_trump20, misreport_trump20)


# Plot 1

# The errorbars overlapped, so use position_dodge to move them horizontally
pd <- position_dodge(0.5) # move them .05 to the left and right


plot_data20$model <- factor(plot_data20$model, levels = c("Validated","Self-reported","Misreported"))


plot2020 <- ggplot(plot_data20, aes(x = x, y = predicted, colour = model)) + 
    geom_point(aes(shape = model), position = pd, size=3.5) + 
    guides(shape=FALSE) + 
    geom_errorbar(
    aes(
      ymin = conf.low, ymax = conf.high,
      colour = model, width = 0.0
    ), size=1, position = pd) + theme_light() + 
    geom_hline(yintercept = 0.1944296, linetype="dashed",
             color = "red", size = 0.8) +
    scale_y_continuous(breaks=seq(0, 0.55, by = 0.1), limits=c(0, 0.55)) + 
  scale_x_continuous(breaks=c(0, 1), 
                     labels = c("0" = "Non Trump Supporter", "1" = "Trump Supporter"), expand = c(0, 0))  +
  xlab("") +
  ylab("Predicted values for Republicans") + 
  scale_shape_manual(values = c(15, 17, 18)) + 
    guides(colour = guide_legend(override.aes = list(shape = c(15, 17, 18)))) +
    theme_bw() + theme(text = element_text(size = 15), 
                     panel.grid.minor = element_blank(), legend.position = "bottom", 
                     axis.title = element_text(size = 12),
                     axis.text.y = element_text(size = 12),
                     axis.text.x = element_text(size = 12),
                     legend.text = element_text(size = 12),
                     legend.title=element_text(size=12)) +
  scale_color_discrete(name="Vote by mail") 

plot2020

# make it black and white 
plot2020 + scale_colour_manual(name=" ", 
                               values = c("Validated" = "#000000", "Self-reported" = "#000000", "Misreported" = "#000000"))


# save graph on overleaf
ggsave(filename = "/Users/Enrijeta/Dropbox/Apps/Overleaf/Leader Motivated Behavior/plots/plot20_rep.pdf", width = 6, height = 6)


##==================##
##==================##
# COVID-19 Exposure ##
##==================##
##==================##

##=====================##
# OLS Regression Models #
##=====================##

m2020_v_covid <- lm(votemethod_validated ~ trump*covid19 + pid_3cat + ideo3 + female + edulevel + 
                                           race_cat + age + married + income + church + news +
                                           noexcuse_vbm + allmail + battleground, data = ces20)
summary(m2020_v_covid)


m2020_s_covid <- lm(votemethod_rv ~ trump*covid19 + pid_3cat + ideo3 + female + edulevel + 
                                    race_cat + age + married + income + church + news + 
                                    noexcuse_vbm + allmail + battleground, data = ces20)
summary(m2020_s_covid)



m2020_covid <- lm(vbm_misreport ~ trump*covid19 + pid_3cat + ideo3 + female + edulevel + 
                                  race_cat + age + married + income + church + news + 
                                  noexcuse_vbm + allmail + battleground, data = ces20)
summary(m2020_covid)



#==============#
# Model Output #
#==============#

# variable names
map <- list("trump" = "Trump supporter",
            "covid19" = "COVID-19",
            "trump:covid19" = "Trump supporterXCOVID-19",
            "pid_3cat2" = "Republican",
            "pid_3cat3" = "Democrat",
            "noexcuse_vbm" = "No excuse absentee",
            "allmail" = "All mail",
            "battleground" = "Battleground",
            "ideo32" = "Conservative",
            "ideo33" = "Liberal",
            "female" = "Female",
            "edulevel" = "Education",
            "race_cat1" = "Black",
            "race_cat2" = "Hispanic",
            "race_cat3" = "Other",
            "age" = "Age",
            "married" = "Married",
            "income" = "Income",
            "church" = "Church attendance",
            "news" = "Political interest",
            "(Intercept)" = "(Intercept)")

# Specify list of models 
models = list(m2020_v_covid, m2020_s_covid, m2020_covid)

# Specify caption 
caption = "Estimating VBM in 2016 and 2020: Validated VBM, Self-reported VBM, and Misreported VBM (COVID-19)"


texreg(models,
       custom.coef.map = map,
       caption = caption,
       caption.above = TRUE,
       digits = 3,
       dcolumn = FALSE,
       stars = numeric(0), 
       override.se = list(
         summary(m2020_v_covid)$coef[,4],
         summary(m2020_s_covid)$coef[,4],
         summary(m2020_covid)$coef[,4]
       )
)


##===============================##
##===============================##
# Appendix: Validated Voters Only #
##===============================##
##===============================##

##===============================##
# Logistic Regression Models 2016 #
##===============================##

m2016_v <- glm(votemethod_validated ~ trump + pid_3cat + ideo3 + female + edulevel + 
                                      race_cat + age + married + income + church + news +
                                      noexcuse_vbm + allmail + battleground, data = ces16, family = binomial())
summary(m2016_v)


m2016_s <- glm(votemethod_rv ~ trump + pid_3cat + ideo3 + female + edulevel + 
                               race_cat + age + married + income + church + news + 
                               noexcuse_vbm + allmail + battleground, data = ces16, family = binomial())
summary(m2016_s)

m2016 <- glm(vbm_misreport ~ trump + pid_3cat + ideo3 + female + edulevel + 
                             race_cat + age + married + income + church + news + 
                             noexcuse_vbm + allmail + battleground, data = ces16, family = binomial())
summary(m2016)



##=========================================================##
# Logistic Regression Models for Vote Method Misreport 2020 #
##=========================================================##

m2020_v <- glm(votemethod_validated ~ trump + pid_3cat + ideo3 + female + edulevel + 
                                      race_cat + age + married + income + church + news + 
                                      noexcuse_vbm + allmail + battleground, data = ces20, family = binomial())
summary(m2020_v)


m2020_s <- glm(votemethod_rv ~ trump + pid_3cat + ideo3 + female + edulevel + 
                               race_cat + age + married + income + church + news + 
                               noexcuse_vbm + allmail + battleground, data = ces20, family = binomial())
summary(m2020_s)



m2020 <- glm(vbm_misreport ~ trump + pid_3cat + ideo3 + female + edulevel + 
                             race_cat + age + married + income + church + news + 
                             noexcuse_vbm + allmail + battleground, data = ces20, family = binomial())
summary(m2020)



##==========================##
# Logistic Regression Models #
##==========================##

all_v <- glm(votemethod_validated ~  trump + pid_3cat + ideo3 + female + edulevel + 
                                     race_cat + age + married + income + church + news + 
                                     noexcuse_vbm + allmail + battleground + year, 
                                     data = ces_all, family = binomial())
summary(all_v)


all_s <- glm(votemethod_rv ~  trump + pid_3cat + ideo3 + female + edulevel + 
                              race_cat + age + married + income + church + news + 
                              noexcuse_vbm + allmail + battleground + year, 
                              data = ces_all, family = binomial())
summary(all_s)


all <- glm(vbm_misreport ~ trump + pid_3cat + ideo3 + female + edulevel + 
                           race_cat + age + married + income + church + news + 
                           noexcuse_vbm + allmail + battleground + year, 
                           data = ces_all, family = binomial())
summary(all)



##=====================##
# DiD Regression Models #
##=====================##

all_vd <- glm(votemethod_validated ~  trump*year + pid_3cat + ideo3 + female + edulevel + 
                                      race_cat + age + married + income + church + news + 
                                      noexcuse_vbm + allmail + battleground, data = ces_all, 
                                      family = binomial())
summary(all_vd)


all_sd <- glm(votemethod_rv ~  trump*year + pid_3cat + ideo3 + female + edulevel + 
                               race_cat + age + married + income + church + news + 
                               noexcuse_vbm + allmail + battleground, data = ces_all, family = binomial())
summary(all_sd)


alld <- glm(vbm_misreport ~ trump*year + pid_3cat + ideo3 + female + edulevel + 
                            race_cat + age + married + income + church + news + 
                            noexcuse_vbm + allmail + battleground, data = ces_all, family = binomial())
summary(alld)



#==============#
# Model Output #
#==============#

# variable names
map <- list("trump" = "Trump supporter",
            "pid_3cat2" = "Republican",
            "pid_3cat3" = "Democrat",
            "noexcuse_vbm" = "No excuse absentee",
            "allmail" = "All mail",
            "battleground" = "Battleground",
            "year" = "Year 2020",
            "trump:year" = "TrumpXYear 2020",
            "ideo32" = "Conservative",
            "ideo33" = "Liberal",
            "female" = "Female",
            "edulevel" = "Education",
            "race_cat1" = "Black",
            "race_cat2" = "Hispanic",
            "race_cat3" = "Other",
            "age" = "Age",
            "married" = "Married",
            "income" = "Income",
            "church" = "Church attendance",
            "news" = "Political interest",
            "(Intercept)" = "(Intercept)")

# Specify list of models 
models = list(m2016_v, m2020_v, all_v, all_vd,
              m2016_s, m2020_s, all_s, all_sd,
              m2016,m2020, all, alld)

# Specify caption 
caption = "Logistic Regression Model for Voting by Mail 2016 and 2020: Validated Voter, Self-reported Validated Voter, and Misreported Validated Method"


texreg(models,
       custom.coef.map = map,
       caption = caption,
       caption.above = TRUE,
       digits = 3,
       dcolumn = FALSE,
       stars = numeric(0),  
       override.se = list(
         summary(m2016_v)$coef[,4],
         summary(m2020_v)$coef[,4],
         summary(all_v)$coef[,4],
         summary(all_vd)$coef[,4],
         summary(m2016_s)$coef[,4],
         summary(m2020_s)$coef[,4],
         summary(all_s)$coef[,4],
         summary(all_sd)$coef[,4],
         summary(m2016)$coef[,4],
         summary(m2020)$coef[,4],
         summary(all)$coef[,4],
         summary(alld)$coef[,4]
       )
)



##===================================================##
##===================================================##
# Appendix: Self-Reported Method, Non-Voters Included #
##===================================================##
##===================================================##

##=====================##
# OLS Regression Models #
##=====================##

m2016_s <- lm(votemethod_r ~ trump + pid_3cat + ideo3 + female + edulevel + 
                             race_cat + age + married + income + church + news + 
                             noexcuse_vbm + allmail + battleground, data = ces16)
summary(m2016_s)


m2020_s <- lm(votemethod_r ~ trump + pid_3cat + ideo3 + female + edulevel + 
                             race_cat + age + married + income + church + news + 
                             noexcuse_vbm + allmail + battleground, data = ces20)
summary(m2020_s)


all_s <- lm(votemethod_r ~  trump + pid_3cat + ideo3 + female + edulevel + 
                            race_cat + age + married + income + church + news + 
                            noexcuse_vbm + allmail + battleground + year, data = ces_all)
summary(all_s)


all_sd <- lm(votemethod_r ~  trump*year + pid_3cat + ideo3 + female + edulevel + 
                             race_cat + age + married + income + church + news + 
                             noexcuse_vbm + allmail + battleground, data = ces_all)
summary(all_sd)



#==============#
# Model Output #
#==============#

# variable names
map <- list("trump" = "Trump supporter",
            "pid_3cat2" = "Republican",
            "pid_3cat3" = "Democrat",
            "noexcuse_vbm" = "No excuse absentee",
            "allmail" = "All mail",
            "battleground" = "Battleground",
            "year" = "Year 2020",
            "trump:year" = "TrumpXYear 2020",
            "ideo32" = "Conservative",
            "ideo33" = "Liberal",
            "female" = "Female",
            "edulevel" = "Education",
            "race_cat1" = "Black",
            "race_cat2" = "Hispanic",
            "race_cat3" = "Other",
            "age" = "Age",
            "married" = "Married",
            "income" = "Income",
            "church" = "Church attendance",
            "news" = "Political interest",
            "(Intercept)" = "(Intercept)")

# Specify list of models 
models = list(m2016_s, m2020_s, all_s, all_sd)

# Specify caption 
caption = "Voting by Mail 2016 and 2020: Self-reported Voter (non-voters included)"


texreg(models,
       custom.coef.map = map,
       caption = caption,
       caption.above = TRUE,
       digits = 3,
       dcolumn = FALSE,
       stars = numeric(0),  
       override.se = list(
         summary(m2016_s)$coef[,4],
         summary(m2020_s)$coef[,4],
         summary(all_s)$coef[,4],
         summary(all_sd)$coef[,4]
       )
)





##=============================================================##
# Appendix: Pooled OLS Regression Models with Year Interactions #
##=============================================================##

all_v <- lm(votemethod_validated ~  trump*year + pid_3cat*year + ideo3*year + female*year + edulevel*year + 
                                    race_cat*year + age*year + married*year + income*year + church*year + 
                                    news*year + noexcuse_vbm*year + allmail*year + battleground*year, data = ces_all)
summary(all_v)



all_s <- lm(votemethod_rv ~  trump*year + pid_3cat*year + ideo3*year + female*year + edulevel*year + 
                             race_cat*year + age*year + married*year + income*year + church*year + 
                             news*year + noexcuse_vbm*year + allmail*year + battleground*year, data = ces_all)
summary(all_s)



all <- lm(vbm_misreport ~ trump*year + pid_3cat*year + ideo3*year + female*year + edulevel*year + 
                          race_cat*year + age*year + married*year + income*year + church*year + news*year + 
                          noexcuse_vbm*year + allmail*year + battleground*year, data = ces_all)
summary(all)



#==============#
# Model Output #
#==============#

# variable names
map <- list("trump" = "Trump supporter",
            "pid_3cat2" = "Republican",
            "pid_3cat3" = "Democrat",
            "noexcuse_vbm" = "No excuse absentee",
            "allmail" = "All mail",
            "battleground" = "Battleground",
            "year" = "Year 2020",
            "trump:year" = "TrumpXYear 2020",
            "ideo32" = "Conservative",
            "ideo33" = "Liberal",
            "female" = "Female",
            "edulevel" = "Education",
            "race_cat1" = "Black",
            "race_cat2" = "Hispanic",
            "race_cat3" = "Other",
            "age" = "Age",
            "married" = "Married",
            "income" = "Income",
            "church" = "Church attendance",
            "news" = "Political interest",
            "year:pid_3cat2" = "RepublicanXYear 2020",
            "year:pid_3cat3" = "DemocratXYear 2020",
            "year:ideo32" = "ConservativeXYear 2020",
            "year:ideo33" = "LiberalXYear 2020",
            "year:female" = "FemaleXYear 2020",
            "year:edulevel" = "EducationXYear 2020",
            "year:race_cat1" = "BlackXYear 2020",
            "year:race_cat2" = "HispanicXYear 2020",
            "year:race_cat3" = "OtherXYear 2020",
            "year:age" = "AgeXYear 2020",
            "year:married" = "MarriedXYear 2020", 
            "year:income" = "IncomeXYear 2020",
            "year:church" = "Church attendanceXYear 2020",
            "year:news" = "Political interestXYear 2020",
            "year:noexcuse_vbm" = "No excuse absenteeXYear 2020",
            "year:allmail" = "All mailXYear 2020",
            "year:battleground" = "BattlegroundXYear 2020",
            "(Intercept)" = "(Intercept)")

# Specify list of models 
models = list(all_v, all_s, all)

# Specify caption 
caption = "Estimating VBM in 2016 and 2020: Validated VBM, Self-reported VBM, and Misreported VBM"


texreg(models,
       custom.coef.map = map,
       caption = caption,
       caption.above = TRUE,
       digits = 3,
       dcolumn = FALSE,
       stars = numeric(0),  
       override.se = list(
         summary(all_v)$coef[,4],
         summary(all_s)$coef[,4],
         summary(all)$coef[,4]
       )
)








